phaseportdat.f <- function (dataset, xv, yv, rangeX, rangeY, f, entidx, plot.limit = 0, ps = F, xlab=xv, ylab=yv, legend=T) 
{
  require(pracma)
  require(MASS)
  procdata <- preprocess_data(2, dataset[,xv], dataset[,yv])
  xwide <- procdata$xwide
  ywide <- procdata$ywide
  tmpx <- rangeX
  tmpy <- rangeY
  xmin = min(tmpx)
  xmax = max(tmpx)
  ymin = min(tmpy)
  ymax = max(tmpy)
  rgrid = 25
  y1 <- linspace(xmin, xmax, rgrid)
  y2 <- linspace(ymin, ymax, rgrid)
  qscale <- range(y1)/range(y2)
  mgrid <- meshgrid(y1, y2)
  x <- mgrid[[1]]
  y <- mgrid[[2]]
  u <- matrix(0, rgrid, rgrid)
  v <- matrix(0, rgrid, rgrid)
  t = 0
  for (i in 1:length(x)) {
    Yprime <- f(t, rbind(x[i], y[i]))
    u[i] <- Yprime[1]
    v[i] <- Yprime[2]
  }
  datasubset <- dataset[!is.na(dataset[,xv])&!is.na(dataset[,yv]),c(xv,yv)]
  kde <- kde2d(datasubset[,1],datasubset[,2],lims=c(rangeX,rangeY))
  u[t(kde$z<plot.limit)] <- NA
  v[t(kde$z<plot.limit)] <- NA
  if (ps) {
    dev.set(1)
    postscript(paste("bdynsys_",xv,"_",yv,".eps",sep=""), width = 5, height = 5, onefile = FALSE, horizontal = FALSE, paper = "special", colormodel = "rgb")#, family = "ComputerModern")
  }
  #layout(cbind(1, 2), widths=c(7, 1))
  #par(mar=c(5.1, 4.1, 4.1, 2.1))
  plot(x, y, col = "white", xlab = xlab, ylab = ylab)
  grid(col = "white")
  quiver(x, y, u, v, scale = 2, angle = 10, length = 0.04, 
         lwd = 1.0, col = "gray")
  rainbow = rainbow(length(entidx))
  for (i in seq_along(entidx)) {
    matplot(xwide[entidx[i], ], ywide[entidx[i], ], type = "l", ldw = 2, col = rainbow[i], add = TRUE)
    points(xwide[entidx[i], 1], ywide[entidx[i], 1], pch = 20, cex = 1.2, col = rainbow[i])
  }
  #par(mar=c(0, 0, 0, 0))
  #plot.new()
  if (legend)
    legend("topleft",xpd=T, inset=c(0,0), bg = "white", cex=.75, pt.cex=.75, legend = entidx, lwd = 1, col = rainbow)
  if (ps)
    dev.off()
}